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1.  Introduction 

The  parameterization  of  turbulence  in  meso-  and  synoptic-scale  numerical  weather  prediction 
(NWP)  models  is  often  accomplished  by  packaging  together  a  large  number  of  highly-specialized 
empirical  approximations  (Anthes,  1983).  Some  of  these  components  handle  the  boundary  layer, 
some  model  clear  air  turbulence  aloft,  and  others  approximate  moist  convection  .  Yet,  these  all 
parameterize  the  same  one  physical  phenomenon  —  turbulence. 

For  example,  the  Penn  State/NCAR  mesoscale  model  (Anthes  et  al.,  1987)  uses  K-theory  vertical 
diffusion,  Blackadar  mixed-layer  parameterization,  bulk  neutral  and  stable  boundary  layer 
parameterizations,  dry  convective  adjustment,  Kuo  cumulus  parameterization,  and  fourth-order 
Smagorinsky  horizontal  diffusion.  Since  many  of  these  schemes  are  applied  to  the  same  grid  points 
within  the  model,  complex  interactions  can  occur  between  the  parameterizations.  To  add  to  the 
confusion,  some  of  these  schemes  are  "tuned"  above  realistic  levels  to  damp  numerical  noise. 

Our  hypothesis  is  that  numerical  weather  forecasts  can  be  improved  when  one  unified  turbulence- 
closure  approximation  is  utilized  for  the  whole  model  domain  at  all  times,  instead  of  using  the 
collection  of  specialized  empirical  approximations.  This  approach  also  requires  that  we  divorce  the 
parameterization  of  the  physics  of  turbulence  from  the  nonphysical  smoothing  necessary  for  numerical 
stability.  In  this  paper,  we  use  transilient  turbulence  theory  (TT)  as  the  unified  turbulence  scheme,  and 
a  6*^-order  implicit  tangent  filter  for  numerical  stability.  The  details  of  TT  are  discussed  in  section  2. 

We  selected  as  a  host  model  the  version  of  the  Penn  State/NCAR  three-dimensional  hydrostatic 
primitive  equation  limited-area  mesoscale  prediction  model  having  fifteen  vertical  sigma  levels  (Anthes 
and  Warner,  1978).  Fig  1  is  a  schematic  showing  that  momentum,  temperature  and  moisture  variables 
are  defined  at  the  half-sigma  level,  while  vertical  velocity  is  defined  on  the  full  sigma  levels  (Anthes  et 
al.,  1987).  In  this  model  the  horizontal  grid  included  61  x  46  points  on  a  Lambert  conformal  mapping 
using  a  horizontal  grid  spacing  of  A)i=Ay=80  km.  The  time  step  is  120  s.  Other  model 
configurations  include  a  nudging  horizontal  boundary  condition  and  a  two  layer  surface  slab 
formulation. 

Although  rarely  stated  in  the  open  literature,  most  nxxielers  will  admit  that  the  various 
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paiameterizations  that  compose  a  mesoscale  model  are  intertwined  in  complex  ways.  Changing  or 
tuning  the  parameters  of  one  parameterization  affect  the  success  of  the  other  parameterizations  in  the 
model  By  changing  to  one  unified  turbulence  thecffy,  we  hope  to  avoid  that  problem.  While  we  were 
untangling  and  removing  the  myriad  small  specialized  turbulence  parameterizations,  we  discovered  an 
impr(^)er  application  of  similarity  theory  fcM*  the  estimation  of  surface  fluxes.  Our  corrected  surface- 
flux  scheme  that  incorporates  a  thin  nonturbulent  viscous  sublayer  is  discussed  in  section  3. 

To  remove  the  numerical  noise  such  as  is  usually  present  in  a  fmite-difference-type  of  mesoscale 
model,  a  6^h_order  implicit  tangent  filter  (Raymond,  1988)  is  installed.  This  filter  is  very  selective  and 
enables  numerical  noise  between  2Ax  and  4Ax  to  be  removed  without  significantly  altering  meaningful 
meteorological  scales.  Thus,  numerical  stability  is  independent  of  the  turbulence  parameterization 
scheme.  Eietails  are  discussed  in  section  4. 

Several  72  h  forecasts  have  been  run  and  verification  made  for  the  OSCAR  IV  (0000  UTC  22  April 
to  0000  UTC  25  April,  1981,  prepared  by  NCAR,  see  Errico  and  Baumhefner,  1987)  and  the  CAPTEX 
(1200  UTC  24  September  to  1200  UTC  27  September,  1983,  prepared  by  NCAR)  data  sets.  The 
former  data  set  represents  a  spring-time  frontal  situation  in  which  cycltmic-frontogenic  activity  is 
located  initially  over  the  Dakotas  and  an  upper-level  trough  exists  just  to  its  west.  This  system 
propagates  eastward  and  intensifies,  but  several  small  scale  waves  are  also  present  (Errico  and 
Baumhefner,  1987).  The  latter  data  set  has  been  used  to  examine  long-range  transport  and  dispersion 
of  pollution  (Kao  and  Yamada,  1988).  During  this  time  period  high  pressure  dominated  most  of  the 
United  States  but  organized  mesoscale  convection  is  jjresent.  Forecast  results  are  presented  in 
Section  5. 

2.  Transilient  turbulence 

Turbulence  is  a  small-scale  phenomenon  for  which  the  horizontal  structure  cannot  be  resolved  in 
large  or  mesoscale  weather  forecast  models,  but  which  affects  the  overall  forecasts.  Convective 
events  originate  from,  and  are  continually  fueled  by,  small  scale  buoyant  plumes  associated  with 
boundary-layer  processes.  The  problem  that  faces  the  modeler  is  how  to  parameterize  these  turbulent 
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motions  within  the  constraints  of  time,  money  and  computer  power,  while  at  the  same  time  resolving 
scales  many  magnitudes  larger. 

The  modelling  of  turbulence  in  the  boundary  layer  and  in  the  atmosphere  above  remains  one  of  the 
more  difficult  problems  in  the  atmospheric  sciences.  Sometimes  turbulence  and  diffusion  have  been 
construed  as  the  same  process.  In  reality  the  processes  involve  different  time  and  space  scales  and  are 
physically  different.  The  first  attempts  to  model  turbulence  (Boussinesq,  1877)  took  the  diffusive 
formulation  for  the  molecular  viscosity  of  air,  enhanced  it,  and  called  it  an  eddy  viscosity.  By 
specifying  the  eddy-diffusion  coefficient  in  terms  of  known  variables,  the  hydrodynamical  equations 
of  motion  can  be  closed  and  solved.  Attempts  to  find  suitable  parameterizations  for  the  eddy  viscosity 
based  on  bulk  or  large  scale  variables  has  occupied  a  predominant  role  in  early  meteorological 
boundary-layer  research.  Following  upon  Boussinesq's  constant  eddy  viscosity  and  Prandtl's  (1925) 
mixing-length  theory,  many  complicated  turbulence  models  have  been  proposed.  In  addition  to  K 
theory's  first-order  closure  scheme  (Louis,  1979)  there  are  higher-order  closure  formulations  (Zeman, 
1981;  Wyngaard,  1982;  Mellor  and  Yamada,  1982;  Andre,  et  al.,  1987)  and  spectral  theories 
(Heisenberg,  1948)  for  turbulence.  Unfortunately,  the  intrinsically  nonlocal  nature  of  atmospheric 
turbulence  makes  it  difficult  to  apply  these  local  approximations  without  violating  the  laws  of 
thermodynamics  (e.g.,  heat  flowing  from  cold  to  hot,  implied  by  negative  K). 

In  this  study  we  turn  our  attention  to  the  newly  developed  transili.;nt  turbulence  parameterization 
(Stull,  1984).  The  nonlocal  transilient  turbulence  parameterization  has  been  compared  with  the  above 
mentioned  local  approaches  by  Stull  (1984, 1986)  and  to  turbulent  adjustment  procedures  by  Stull  and 
Hasegawa  (1984).  In  many  one  dimensional  tests  (Stull  and  Driedonks,  1987;  Stull  and  Kraus,  1987) 
this  new  formulation  has  yielded  outstanding  results. 

a.  The  nonlocal  nature  of  turbulence 

Turbulence  in  the  atmosphere  is  not  like  diffusion,  but  is  more  like  advection.  Blobs  or  parcels  of 
air  are  bodily  moved  across  finite  distances  by  eddies  of  varying  sizes,  before  the  parcels  (or  portions 
of  them)  mix  with  their  surroundings.  This  picture  is  particularly  valid  for  both  boundary-layer 
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convection  and  convective  clouds  within  the  troposphere.  For  statically  stable  air,  the  same  picture  is 
valid,  although  the  dominant  eddies  are  usually  smaller.  Thus,  a  desirable  feature  of  turbulence 
parameterizations  for  numerical  weather  forecast  models  is  that  mixing  across  all  relevant  scales  be 
adequately  described. 

The  concept  of  mixing  across  finite  distances  is  embodied  in  a  nonlocal  first-order  closure  scheme 
called  transilient  turbulence  theory  (TT).  The  word  "transilient"  is  based  on  a  latin  word  meaning 
"jump  over",  to  convey  the  advective  rather  than  diffusive  nature  of  TT.  For  a  1-D  column  of  grid 
points,  there  is  a  matrix  [cy]  of  mixing  coefficients  (transilient  coefficients)  that  describe  how  much  air 
arriving  at  grid  box  i  came  from  box  j  during  timestep  At.  Grid  points  i  and  j  can  represent  immediate 
neighbors,  or  they  can  be  points  separated  in  space  (see  Fig  2).  Thus,  the  full  range  of  resolvable 
scales  can  be  represented.  In  Fig  2b,  the  curved  arrow  represents  a  mixing  process,  which  is  the  net 
effect  of  a  superposidon  of  actual  eddies  that  contribute  to  mixing  between  the  two  grid  points  during 
timestep  At.  Similar  representations  can  be  made  in  the  horizontal. 

Measurements  of  the  transilient  matrix  have  been  made  by  Ebert  et  al.  (1989)  for  the  case  of  free 
convection  in  the  boundary  layer,  using  equally-spaced  grid  points.  They  ran  a  large-eddy  simulation 
model  and  injected  tracers  to  track  the  sources  and  destinations  of  all  the  air  within  the  model.  It  was 
not  surprising  that  they  found  all  scales  of  motion  to  be  important  to  vertical  mixing,  not  just  the  small- 
scales  modeled  by  traditional  K-theory.  In  addition,  they  found  that  the  matrix  is  asymmetric  for 
unstable  boundary  layers.  The  physical  interpretation  of  any  transilient  matrix  is  discussed  in  detail  by 
them. 

b.  Implementation  of  turbulent  mixing 

Implementation  of  this  scheme  into  a  1-D  forecast  model  (or  into  any  column  of  a  3-D  model) 
requires  parameterization  of  the  transilient  matrix.  The  parameterization  should  be  able  to  respond  to 
changes  of  the  mean  state  caused  by  the  various  body  fOTcings  (radiation,  condensation,  Coriolis 
force,  etc.)  and  boundary  conditions  (surface  heat  flux,  drag,  etc.). 

The  earliest  parameterizations  (Berkowicz  and  Prahm,  1979)  assumed  a  fixed  spectrum  of 
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turbulence  that  was  not  particularly  responsive  to  changes  in  the  mean  state.  Later,  a  responsive 
parameterization  was  developed  by  Stull  (1984)  based  on  a  nonlocal  Richardson  number,  but  it  did  not 
perform  well  for  arbitrary  situations  such  as  radiative  cooling  at  cloud  top  within  the  model  domain. 
More  recently,  Stull  and  Driedonks  (1987)  developed  a  responsive  parameterization  based  on  a 
nonlocal  approximation  to  the  turbulence  kinetic  equation  that  gave  realistic  forecasts  for  unstable, 
neutral,  stable,  and  cloudy  boundary  layers.  Although  this  scheme  yields  only  symmetric  transilient 
matrices,  it  is  very  robust  and  automatically  satisfies  conservation  of  mass  and  state.  However,  the 
parameterized  turbulence  is  biased  to  generate  the  most  intense  mixing  at  the  lower  grid  points,  which 
can  be  corrected  as  is  discussed  later. 

The  responsive  parameterizations  require  that  there  be  some  instability  to  respond  to.  In  addition,  a 
3-D  model  requires  both  horizontal  and  vertical  mixing.  Thus,  each  model  forecast  timestep  must  be 
split  into  five  parts;  (i)  application  to  all  grid  points  of  all  (nonturbulent)  boundary  conditions  and  body 
forcings  including  dynamics  and  thermodynamics  that  could  alter  or  destabilize  the  flow;  (ii) 
computation  of  a  transilient  matrix  for  vertical  mixing  for  each  column  of  grid  points  based  on  this 
destabilized  mean  state;  (iii)  vertical  mixing  of  all  state  variables  in  each  column  based  on  the 
computed  transilient  matrix;  (iv)  computation  of  horizontal  mixing  potentials  for  each  column;  and  (v) 
horizontal  mixing  of  all  state  variables.  The  procedure  is  first  do  step  (i)  over  the  whole  3-D  domain, 
next  do  steps  (ii)  and  (iii)  for  each  column,  and  finally  do  steps  (iv)  and  (v). 

Step  (i)  is  conducted  using  the  existing  physics  of  the  host  model,  with  the  exception  that  all 
turbulence-related  effects  are  excluded.  In  the  Penn  State/NCAR  model,  we  removed  vertical  K-theory 
diffusion,  Smagorinsky's  horizontal  diffusion,  dry  convective  adjustment,  Kuo  subgrid  cumulus 
parameterization,  and  all  the  boundary  layer  paranwterizations.  Note  that  the  flux  of  heat,  moisture, 
and  momentum  between  the  ground  and  the  air  is  not  a  turbulent  process  (unless  the  ground  begins  to 
dance),  and  hence  that  parameterization  is  retained  (see  Section  3). 

The  net  effect  of  step  (i)  is  to  split  each  timestep  into  (at  least)  two  parts:  the  nonturbulent-forcing 
part  that  incorporates  all  the  other  physics,  followed  by  the  turbulence  part.  At  the  beginning  of  any 
timestep,  let  the  state  of  air  within  any  one  column  of  grid  points  be  Sj(t),  where  Sj  represents 
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moisture,  temperature,  or  momentum-component  values  at  each  grid  point  i  within  the  column.  After 
applying  the  physics  of  step  (i),  define  the  altered  state  as  Sj'(t). 
c.  Vertical  mixing 

For  step  (ii),  we  start  with  the  1-D  approach  of  Stull  and  Driedonks  (1987)  to  determine  the 
transilient  matrix  for  any  one  column  of  grid  points.  The  wind  speed  differences  [(AU)ij  and  (AV)ijl 
and  virtual  potential  temperature  difference  [(A0v)ij  1  between  each  pair  (i  j)  of  grid  points  within  the 
column  determines  a  mixing  potential  (Yjj)  between  those  points: 
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for  i>tj,  g  is  gravitational  acceleration,  and  where  (A0v)ij  is  negative  if  the  lower  grid  point  is  wanner 
than  the  higher  point.  The  distance  between  pairs  of  g  id  points  (i.e.,  between  the  centers  of  the  grid 
boxes)  is  Azjj  as  sketched  in  Fig  2b.  The  above  equation  is  also  valid  for  unequally- spaced  grid 
points.  Note  that  a  typographical  error  in  (6)  of  Stull  and  Driedonks  (1987)  resulted  in  the  omission 
of  the  Azij  factor  in  the  last  term  in  square  brackets,  but  that  their  calculations  and  results  were  based 
on  the  correct  formula. 

The  Penn  State/NCAR  mesoscale  model  uses  a  staggered  grid  (Arakawa  B).  Consequently, 
velocities  are  stored  at  locations  horizontally  displaced  from  temperature  and  humidity.  Interpolations 
of  velocity  are  made  to  coincide  to  temperature  locations  before  the  Ajj  differences  in  (1)  are  computed. 
Although  this  procedure  is  required  for  physical  consistency  with  this  particular  host  model,  it 
unfortunately  causes  extra  smoothing  which  can  reduce  the  Ajj  differences  and  alter  the  mixing 
potential. 

Parameter  values  suggested  by  Stull  and  Driedonks  are  also  used  here:  timescale  Tq  =  1000  s, 
critical  Richardson  number  Rc  =  0.21 ,  and  dissipation  parameter  D  =  1 .  Along  any  row  of  the 
mixing-potential  matrix,  values  closer  to  the  main  diagonal  are  increased  if  necessary  to  be  no  smaller 
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than  elements  further  from  the  diagonal  to  eliminate  unrealistic  convective  overturning.  The  internal 
subgrid  scale  mixing  potential,  Yu,  is  then  set  equal  to  the  largest  element  on  its  same  row,  plus  a 
reference  potential  Yref  =  1000.  The  net  result  is  that  we  have  used  4  parameters  (To  ,  Rc,  D,  and  Yj-^f) 
to  parameterize  a  mixing-potential  matrix  with  15x15  unknowns. 

At  this  point,  the  unequally-vertically-spaced  grid  points  in  the  host  model  require  a  modification  of 
the  original  equations  suggested  by  Stull  and  Driedonks.  A  row-norm  (RNj)  for  each  row  in  the 
transilient  matrix  is  defined  by  the  weighted  sum  of  the  mixing  potentials; 


m.  Y. 


j  'J 


(2) 


where  mj  is  the  mass  of  air  within  grid  layer  j.  Alternately,  one  can  use  the  relative  mass  (or  Aa 
within  a  sigma-coordinate  model  such  as  the  Penn  State/NCAR  model)  in  place  of  mj.  The  so-called 
Leo  norm  (i.e.,  the  value  of  the  one  maximum  row)  is  defined  as  a  matrix-norm  IIYIleo  ,  as  in  Stull  and 
Driedonks; 

ilYII  =  max  (RN.)  (3) 

~  i  ' 

The  off-diagonal  transilient  coefficients  are  parameterized  by; 


c..  = 


m. 

J 


Y 

u 


IlYII 


(4a) 


and  the  elements  on  the  main  diagonal  are  given  by; 


c.. 


U 


j  =  i 
j  i 


(4b) 


It  is  important  to  note  that  Cij  ^  Cji,  even  if  Yij  =  Yji,  because  mj  ^  mj.  Thus,  Yjj  is  symmetric. 


7 


while  Cjj  is  asymmetric.  Furthermore,  each  row  of  the  transilient  matrix  sums  to  one,  but  each  column 
does  not.  In  spite  of  this,  total  air  mass  and  state  variables  arc  still  conserved  during  the  mixing 
process,  and  the  mixing  is  assumed  to  represent  an  equal  (but  not  necessarily  total)  exchange  of  mass 
between  the  grid  points  i  and  j.  The  above  scheme  works  well  with  any  arbitrary  grid  spacing, 
including  the  sigma  coordinates  of  the  host  model.  A  more  detailed  explanation,  derivations  of  the 
above  equations,  and  examples  of  transilient  matrices  for  unequally-spaced  grid  points  are  given  in 
Appendix  A. 

Equations  (2)-(4)  above  are  exact  expressions,  based  on  conservation  of  mass;  however,  (1)  is  a 
parameterization  that  is  only  an  approximation  to  nature.  Recently,  Stull  (1989)  compared  the 
parameterized  matrices  with  the  matrices  observed  by  Ebert  et  al.  (1989),  and  found  a  systematic  bias 
in  (1).  The  parameterized  matrices  tend  to  cause  erroneously  large  amounts  of  mixing  at  the  bottom  of 
each  turbulent  domain,  compared  to  the  mixing  between  grid  points  elsewhere  within  the  domain,  even 
for  equally-spaced  grid  points.  A  similar  conclusion  was  discovered  earlier  by  Chrobok  (1988)  based 
on  boundary-layer  forecasts  made  with  (1).  For  this  reason,  we  designed  a  weighting  scheme  to 
counteract  the  unrealistic  height-dependent  mixing.  For  simplicity,  we  divided  (4a)  by  mj,  which 
works  for  the  Penn  State/NCAR  model  only  because  the  grid  layer  thickness  increases  monotonically 
with  height.  The  net  effect  is  to  reduce  the  mixing  into  lower-altitude  grid  layers  while  increasing  it  in 
the  upper  layers.  Although  this  correction  scheme  was  used  here  as  a  quick  fix,  a  better  approach 
would  have  been  to  refine  (1)  into  a  more  physical ly-realistic  parameterization. 

During  any  one  timestep,  the  operations  of  equations  (1)  -  (4)  described  above  can  be  applied  to 
one  column  at  a  time,  until  the  whole  domain  of  the  host  model  has  been  turbulently  adjusted. 
Occassionally  we  have  stored  the  transilient  matrix  for  each  column  for  post  analysis  of  turbulence 
characteristics;  however,  storage  limitations  on  the  computer  precluded  extensive  use  of  this  approach. 
Instead,  we  usually  stored  only  ye.s/no  information  on  the  turbulence  state  of  each  grid  point.  These 
vertical  cross  sections  turbulence  state  are  useful  to  highlight  the  boundary  layer,  frontal  zones,  and 
clear-air  turbulence  patches.  Horizontal  cross-sections  show  unstable  air  masses  and  frontal  zones. 
Also,  this  information  of  turbulent  subdomains  is  important  input  to  the  horizontal  mixing  scheme 
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described  in  the  next  subsection.  As  a  compromise  between  the  storage  of  full  turbulence  details  and 
limited  computer  storage,  we  sometimes  save  only  the  Cjj  value  for  each  grid  point  i,  because  it  not 
only  indicates  if  the  point  was  turbulent  (Cjj  <  1.0),  but  the  amount  of  mixing  is  indicated  by  the 
magnitude  of  1.0-  cjj. 

Next,  step  (iii)  is  accomplished  with  simple  matrix  multiplication  applied  to  each  of  the  mean  field 
state  variables  in  a  column: 


S.  (t+At)  =  ^  c..  (t.  At)  S/(t )  (5) 

j  =  1 

Steps  (ii)  and  (iii)  are  applied  together,  one  column  at  a  time. 
d.  Horizontal  mixing 

For  step  (iv),  we  explicitly  assume  that  turbulence  is  a  three-dimensional  phenomenon.  If  there  is 
(resolved)  turbulent  muting  in  the  vertical,  then  there  is  also  mixing  by  the  same  (unresolved)  eddies  in 
the  horizontal.  If  no  turbulence  has  been  found  in  a  portion  of  a  column  based  on  the 
parameterizations  of  the  previous  subsection,  then  we  assume  that  horizontal  mixing  cannot  occur  their 
either.  For  each  height,  we  assume  that  horizontal  mixing  between  two  neighboring  columns  occurs 
only  if  both  columns  have  turbulence  at  that  height  (see  Fig  3).  Also,  we  are  not  concerned  with  non¬ 
neighboring  columns,  because  the  column  width  is  so  large  (80  km)  and  the  timestep  of  the  host  model 
is  so  small  (120  s)  that  turbulence  can  not  physically  mix  beyond  its  neighbor  during  one  timestep. 
This  horizontal  mixing  between  neighbors  at  first  glance  looks  similar  to  K-theory,  except  that  here  we 
allow  horizontal  mixing  only  between  neighbors  that  are  turbulent. 

The  procedure  is  to  examine  each  column  one  at  a  time.  Define  the  column  of  interest  by  sub.script 
i  =  0,  and  each  of  its  (up  to)  four  neighbors  on  the  same  sigma  surface  by  sub.script  i  =  1  to  4.  For 
each  column,  we  then  examine  each  height  one  at  at  time.  Based  on  the  existence  of  turbulence  and  the 
vertical  size  of  the  contiguous  turbulent  domain,  we  define  horizontal  mixing  potentials  Vyj  at  that  one 
height. 
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We  assume  the  hcnizontal  mixing  potential  is  equal  to  the  fraction  of  a  grid  box  width  into  which 
uirbulence  can  transport  state  variables  from  a  neighboring  box.  Using  simple  eddy-size  arguments, 
this  fraction  is  propcational  to  the  lateral  distance  (h^o)  that  turbulence  can  mix  during  one  eddy-cycle, 
times  the  number  of  eddy-cycles  (At/t)  that  occur  during  one  timestep,  divided  by  the  horizontal  size  of 
the  grid  box.  Ax: 

h  .  At 

V  _  ^ 


In  this  parametehzation,  we  have  assumed  for  computational  efficiency  that  turbulent  eddies  are 
isotropic,  and  that  the  lateral  mixing  distance  is  proportional  to  the  size  of  the  largest  eddies,  which  in 
turn  is  the  vertical  size  (h^o)  of  the  contiguous  regions  of  turbulence  that  are  simultaneously  present  in 
the  neighboring  columns.  The  effective  timescale  for  one  eddy  cycle,  x,  is  unknown,  and  is  treated  as 
a  parameter  to  be  tuned  during  trial  forecasts.  We  found  x  =  100  s  to  work  best,  given  the  NCAR/ 
PSU  timestep  of  120  s.  The  subgrid-scale  mixing  potential  (Vqo)  is  defined  by  (6),  with  hjo  replaced 
by  the  total  vertical  domain  of  the  host  model.  Furthermore,  Voi  =  0  at  any  height  if  either  the  center 
column  (o)  or  its  neighboring  column  (i)  is  not  turbulent  at  that  height. 

Step  (v)  can  now  be  performed  at  each  grid  point  within  each  column  by  applying  the  Lj  norm 
(shown  in  the  denominator  of  (7)]: 


So(after  horizontal  mixing) 


i 


i  =  n 


S.  (before  horiz.  mixing) 


i  =  0 


(7) 


The  results  of  horizontal  mixing  are  stored  in  a  separate  array  until  the  computations  have  been 
completed  for  the  whole  domain.  Thus  the  mixing  potentials  at  a  specified  grid  point,  as  determined 
by  applying  (6)  and  (7),  are  based  on  a  5-point  stencil  (see  Fig  4).  As  a  side  note,  (7)  can  theoretically 


be  put  in  matrix  form,  and  combined  with  (5)  to  yield  one  equation  that  can  be  inverted  to  yield  a 
different  c  matrix  that  accomplishes  both  the  vertical  and  horizontal  mixing. 

Some  initial  forecast  tests  were  made  with  no  horizontal  mixing.  In  general,  the  solutions  were 
numerically  stable  (both  with  and  without  the  filter  described  in  section  4);  however,  frontal 
boundaries  became  too  sharp  as  the  forecast  progressed.  When  horizontal  mixing  was  added,  the 
frontal  zones  became  more  realistic.  Also,  horizontal  mixing  alters  the  temperature  and  wind  profiles 
in  some  of  the  columns,  which  in  turn  alters  the  dynamic  stability  and  vertical  mixing.  The  combined 
nonlocal  vertical  and  horizontal  mixing  are  very  strongly  interdependent.  The  result  is  a  quasi-3-D 
turbulence  scheme  that  automatically  forecasts  boundary-layer  evolution  and  patchy  clear  air 
turbulence. 

e.  Computational  efficiency 

We  now  examine  the  number  of  calculations  involved  in  the  turbulence  parameterization  during  any 
one  time  step  for  a  grid  with  N  x  M  horizontal  and  L  vertical  discrete  points. 

The  staggering  of  the  grid  in  the  Penn  State/NCAR  model  (Arakawa  B  grid)  requires  that  some 
interpolation  be  done  in  the  calculation  of  the  mixing  or  transilient  coefficients  in  both  the  vertical  and 
horizontal  turbulence  schemes.  Specifically,  the  horizontal  velocities  are  linearly  interpolated  to 
temperature  grid  point  locations.  These  interpolated  values  are  then  used  in  the  calculation  of  the 
mixing  potential  for  the  temperature  and  mixing  ratio.  The  transilient  coefficients  associated  with  the 
temperature  and  mixing  ratio  fields  are  then  retained  and  linearly  interpolated  to  the  staggered  grid 
locations  to  give  the  transilient  coefficients  associated  with  the  horizontal  wind.  Horizontal  boundaries 
are  assigned  no-mixing  potentials.  Because  of  the  staggered  grid  many  extra  calculations  are  required. 

Bookkeeping  processes  require  147  NML  arithmetic  operations;  i.e.,  operations  requiring  addition, 
subtraction,  division  or  multiplication  of  a  real  or  integer  quantity.  Determination  of  the  mixing 
potentials  for  use  in  (5)  can  require  between  189  and  250  NML  operations  depending  on  the  amount 
of  turbulence.  Matrix  multiplication  requires  120  NML  evaluations.  All  total,  the  greatest  possible 
number  of  calculations  required  for  an  application  of  (5)  to  the  entire  domain  is  517  NML  operations 
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not  counting  a  maximum  of  9  NML  determination  of  a  maximum  value.  If  every  point  is  turbulent 
the  horizontal  mixing  formulation  requires  179  NML  operations.  This  give  somewhere  in  the 
neighborhood  of  700  NML,  or  less,  arithmetic  operations  fw  the  combined  vertical  and  horizontal 
turbulence  calculations.  This  is  a  substantial  number  of  arithmetic  operations,  requiring  more  than  a 
three-fold  increase  in  the  total  computing  time  to  make  a  72  hr  forecast,  provided  the  turbulence 
paranseterization  is  applied  at  every  time  step. 

We  will  show  in  Section  5  that  the  transilient  turbulence  parameterization  can  increase  the 
magnitude  of  the  omega  field  by  a  factor  of  two,  and  therefore  models  severe-storm  precipitation 
events  better.  It  is  well  known  that  the  magnitude  of  the  vertical  velocity  is  dependent  on  the  grid 
resolution.  Thus,  we  would  expected  similar  behavior  without  TT  if  the  grid  resolution  was  enhanced 
by  a  factor  of  two  or  more.  If  this  was  done,  small  scale  horizontal  structures  would  be  better 

represented  and  the  omega  field  would  benefit  since  the  horizontal  divergence  is  scale  dependent. 
Doubling  of  the  grid  resolution,  however,  would  cause  an  eight-fold  increase  in  computing  time. 

When  measured  against  this  increase,  the  transilient  3-foId  increase  provides  an  economic  alternative. 

A  significant  reduction  (about  one-fifth)  in  the  total  number  of  required  calculations  could  be 
obtained  if  the  grid  was  not  staggered.  Likewise  improved  vcctorization  would  help  but  the 
parameterization  is  strongly  centered  around  the  vertical  dimension,  which  being  the  smallest,  limits 
the  potential  for  improvement.  We  are  optimistic  that  future  research  will  find  ways  to  significantly 
reduce  the  computational  expense.  However,  the  final  authority  that  determines  whether  the 
parameterization  is  worth  the  expense  is  the  quality  of  the  forecasts. 

3.  Surface  heat  flux  with  molecular  effects  included 

One  of  the  major  roles  of  turbulence  is  to  transport  heat  and  moisture  from  the  surface  to  the  rest  of 
the  atmosphere.  Turbulence  cannot  accomplish  its  task  until  molecular  processes  have  first  transponed 
the  heat  and  moisture  from  the  nonturbulent  surface  into  the  lowest  few  millimeters  of  air.  Thus, 
transport  from  the  ground  to  the  air  is  controlled  by  both  molecular  and  surface-layer  physics. 

Flux-profile  relationships  for  the  surface  layer  empirically  describe  only  the  turbulent  contribution 
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to  flux.  Given  measurements  of  potential  temperature  difference  between  two  heights  within  the 
turbulent  surface  layer,  for  example,  a  heat  flux  can  be  estimated.  It  is  physically  inappropriate  to  use 
a  temperature  difference  between  the  ground/sea  surface  and  the  air  in  such  a  surface-layer 
parameterization.  We  show  below  how  surface  layer  and  molecular  layer  parameterizations  can  be 
combined  to  estimate  flux  from  the  ground-to-air  temperature  difference. 

The  heat  flux  used  in  the  surface-layer  turbulence  calculations  within  the  Penn  State/NCAR  model 
is  computed  from  similarity  theory  (Businger  et  al.,  1971): 

H3  =  -pCpku.T,  .  (8) 

where  is  given  by 

T,  =  (0-e,)/(ln(z^Zo)-VH).  (9) 

The  parameter  p  is  the  density,  k  is  the  von  Karman  constant  of  0.4,  Cp  is  the  specific  heat  of  air  at 
constant  pressure,  9^  is  the  potential  temperature  at  a  nominal  "surface"  height  z^  while  0^  is  the  air 
potential  temperature  at  the  lowest  model  level  z,.  The  friction  velocity,  u,,  is  expressed  as 

u,=Max[  u  ,  k  V  /  (In(Za/Zo)  -  w  )  ] ,  (10) 

®  M 

where  u,^  =  0.1  ms'^  and  V=(Vg^  +  Here  Vg  is  the  wind  speed  at  height  Zg  and  V^,  is  a 

convective  velocity  (Anthes  et  al.,  1987).  The  toughness  height  is  z^.  The  nondimensional  stability 
parameters  and  are  a  function  of  the  bulk  Richardson  number  (Anthes  et  al.,  1987). 
Combining  (8)  thru  (10)  yields 


H,  --Ci(0g-0J. 


(11) 
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The  original  parameterization  erroneously  uses  the  ground  skin  or  sea-surface  temperature  for  85, 
which  violates  the  assumptions  of  the  similarity  theory.  Instead,  we  will  assume  that  6^  is  the 
temperature  in  the  air  at  the  bottom  of  the  turbulent  surface  layer,  which  corresponds  with  the  height 
of  the  top  of  the  nonturbulent  viscous  microlayer  (Stull,  1988). 

In  the  molecular  layer,  the  molecular  heat  flux  at  is  given  by  the  diffusion  equation: 

Hsu  =  -P  KH^  (Ts  -  Tg)  /  Zj,  -  -C2(0,  -  ep.  (12) 

Here  Khjx  is  the  molecular  conductivity  in  air,  and  9g  is  the  surface-skin  temperature. 

By  matching  the  heat  fluxes  (H  s  Hj  =  )  and  65  between  the  surface  layer  and  microlayer,  we 

can  combine  (1 1)  and  (12)  to  yield: 


H  =  -ci(0,-0g)(l+ci/c2)-'  (13a) 

or 

H  =  -CiC3(0.-0g)  (13b) 

where 

C3  a  (1+  Cj/C2)-’  .  (14a) 

and 

C1/C2  =  Cp  k^  V  z^  /[KHpln(Za/zo  -  \|r^)ln(z^Zo  -  Vj,)]  (14b) 

We  see  that  (13b)  is  similar  to  (1 1),  except  that  the  additional  factor  C3  appears  when  the  ground 
skin  temperature  is  used  in  place  of  a  surface  air  temperature.  For  typical  values  of  V,  z^,  z^.  and 
(e.g.,  V  =  5  ms-i,  z,  =  50  m,  z^  =  0.01  m  and  between  0.01  and  0.5),  we  obtain  for  C3  a  range  of 
values  between  0.20  and  0.06.  Consequently  it  is  clear  that  the  molecular  layer  reduces  the  heat  flux  in 
(13b)  by  at  least  80%  for  any  given  temperature  difference. 
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Over  land  surfaces,  however,  the  ground  skin  temperature  is  not  fixed,  but  responds  to  the  surface 
energy  balance.  Thus,  a  fixed  rate  of  insolation  will  cause  the  ground  skin  temperature  to  become 
warmer  than  in  the  original  Penn  State/NCAR  model,  resulting  in  a  new  equilibrium  heat  flux  that  is 
nearly  the  same  magnitude  as  in  the  original  model  (i.e.,  it  will  be  5  - 10%  less,  rather  than  80%  less). 
Although  there  is  little  net  change  in  heat  flux,  the  ground  temperature  is  now  significantly 
warmer/cooler  than  the  air  temperature  during  day/night,  and  there  is  a  strong  temperature  gradient 
across  the  microlayer.  In  our  model  runs  we  have  tested  forecasts  where  the  value  of  C3  is  fixed  at 
0.2,  0.1  or  allowed  to  vary  over  a  fixed  range  (0.1  <  C3  <  0.2).  For  all  cases  we  found  (13b)  to  be  a 
more  realistic  and  consistent  parameterization  for  the  heat  flux  at  the  surface  than  the  original 
parameterization. 

The  original  surface  moisture  flux,  E,  parameterization  in  the  Penn  State/NCAR  model  was 
similarly  modified  to  include  the  C3  parameter  associated  with  the  viscous  microlayer  over  land: 

E  =  -CjCj  -^(q,  -  Qg)  (15) 

where  Mg  is  the  moisture  availability  parameter,  q^  is  the  specific  humidity  at  the  lowest  grid  point  in 
the  air,  and  the  ground  moisture  qg  is  taken  as  the  saturation  specific  humidity  at  the  ground  skin 
temperature. 

4.  The  Implicit  Tangent  Filter 

Numerical  forecast  models  normally  require  dissipative  finite  difference  schemes  or  artificially 
enhanced  viscosity  to  maintain  numerical  stability.  The  Penn  State/NCAR  regional  model  is  known  to 
become  unstable  when  the  horizontal  diffusion  is  reduced  below  a  critical  value  (Errico  and 
Baumhefner,  1987).  In  contrast,  we  find  the  TT  version  of  the  model  to  be  numerically  stable,  for  the 
cases  considered,  when  we  removed  all  K-  theory  diffusion.  This  is  true  with  or  without  the 
horizontal  mixing  component  in  the  turbulence  parameterization.  Stull  (1986)  had  proved  that  TT  by 
itself  is  numerically  stable;  however,  we  were  surprised  to  learn  that  the  host  model  remained  stable,  in 
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spite  oi  the  patchy  spOTadic  turbulence  produced  by  TT. 

Nevertheless,  the  sporadic  nature  of  the  forecasted  upper-level  turbulence,  and  the  growth  and 
decay  in  the  boundary  layer  did  introduce  light  to  moderate  amounts  of  numerical  noise  in  the 
horizontal  directions.  Even  though  the  magnitudes  of  these  features  remained  fmite  in  time  it  is  still 
desirable  to  remove  all  scales  that  cannot  be  represented  accurately  by  the  finite  difference  scheme. 
Filters  are  commonly  used  in  numerical  models  to  remove  high  frequency  noise  that  can  not  be 
resolved.  Among  the  most  familiar  are  the  Shapiro  (1970,1975)  and  Shuman  (1957)  filters.  In  this 
study  we  apply  the  sixth  order  low-pass  implicit  tangent  filter  (Raynwnd,  1988).  This  filter  possesses 
some  very  desirable  characteristics,  viz.,  it  has  an  extremely  sharp  roll-off  in  the  amplitude  response 
function,  can  be  applied  near  horizontal  boundaries  and  it  contains  a  filter  parameter  e  (equation  in 
APPENDIX  B)  which  can  be  adjusted  to  give  the  desired  smoothing. 

The  implicit  calculations  require  the  inversion  of  a  banded  diagonally  dominant  matrix.  For  t»der 
two,  the  filter  studied  in  Pepper,  et  al.,  (1979)  is  recovered.  As  the  wder  is  increased  the  roll-off  in  the 
amplitude  response  is  greatly  sharpened.  The  characteristics  of  the  sixth-order  low-pass  implicit 
tangent  filter  are  best  illustrated  by  examining  the  amplitude  response  function  given  in  Appendix  B. 
This  response  is  similar  to  that  given  by  the  recursive  tangent  filters  described  in  Otnes  and  Enochson 
(1978).  However  the  implicit  tangent  filter  is  nevertheless  much  easier  to  woric  with  because  the 
coefficient  weights  are  known  and  boundary  conditions  encountered  in  limited  area  modelling  are 
easily  handled  in  the  implicit  formulation. 

In  Fig.  5  the  response  of  the  filter  is  shown  after  24  hrs  (720  applications  with  e=0.0075)  and  is 
compared  with  the  4th  order  K  theory  horizontal  smoothing  normally  used  in  the  Penn  State/NCAR 
regional  model.  Note  for  wave  number  6Ax  that  the  responses  are  much  closer  to  the  unfiltered  value 
of  1  when  using  the  6th  order  implicit  filter,  whereas  the  4th  order  K  theory  operator  removes  almost 
all  of  these  features,  especially  when  the  maximum  K  value  is  used.  The  filter,  with  e=.0075,  is 
applied  at  every  time  step  in  the  model  to  the  horizontal  wind  velocity  components,  the  temperature  and 
the  mixing  ratio  fields.  The  filter  parameter  e  was  selected  to  give  the  lowest  acceptable  smoothing 
while  still  maintaining  reasonable  smooth  fields.  The  application  of  the  filter  takes  place  after 
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completion  of  the  new  time  step  and  following  the  turbulence  parameterization. 

5.  Discussion  of  Results 

a.  The  OSCAR  TV  case 

The  percent  of  the  total  number  of  horizontal  grid  points  that  contain  turbulent  exchanges,  for  the 
OSCAR  IV  case,  is  illustrated  in  Fig.  6  as  a  function  of  pressure  and  time,  as  predicted  from  a  72  hr 
forecast  with  the  TT  panuneterization.  Note  that  above  the  500  mb  level  approximately  one  percent  or 
less  of  the  total  number  of  the  horizontal  grid  points  are  turbulent.  The  maximum  at  high  levels  occurs 
during  the  strongest  cyclonic  intensification  which  takes  place  between  hours  18  through  48.  In  the 
lowest  levels  of  the  boundary  layer,  more  than  75%  of  the  total  horizontal  grid  points  are  turbulent 
during  the  peak  heating  in  the  diurnal  cycle.  Qearly  there  is  a  significant  diurnal  variation  within  the 
boundary  layer. 

The  horizontal  distribution  of  turbulence  in  the  boundary  layer  is  shown  in  Figs.  7a,b.  In  Fig.  7a 
the  distribution  of  the  turbulence  21  hours  into  the  forecast  (near  3  pm  local  time)  at  a=0.94 
corresponds  closely  with  the  heated  land  mass.  Only  in  the  midwest  (associated  with  the  major 
cyclonic  activity)  is  there  a  turbulence-free  zone.  This  most  likely  occurs  because  of  the  stabilizing 
radiative  and  evaporative  properties  associated  with  the  precipitating  clouds.  Also  note  the  lack  of 
turbulence  over  the  water  surface.  At  higher  levels  in  the  atmosphere  (a=.74)  the  turbulence  is 
confined  to  the  mountainous  regions  as  shown  in  Fig.  7b  . 

Vertical  cross  sections  through  the  center  of  our  region  (taken  west  to  east)  are  presented  in  Figs. 
8a,b  showing  the  potential  temperature  0  (K)  and  the  mixing  ratio 

(g  kg"  1).  Regions  with  turbulence  are  within  or  under  the  wide  solid  line  found  in  the  0  field  in  Fig. 
8a.  The  turbulent  (mixed)  boundary  layer  appears  deepest  over  the  mountains  and  just  ahead  of  the 
front  as  shown  between  grid  points  30  and  45,  which  is  very  realistic.  Some  nirbulence  is  also 
occurring  at  mid-levels  within  the  frontal  zone  above  grid  points  29  and  30.  Note  that  some  folding  in 
the  0  contours  occur  at  mid  levels  between  grid  points  25  and  35  indicating  the  presence  of  the  cold 
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front.  Clearly  the  largest  amount  of  moisture  is  found  in  the  warm  sector.  This  is  reflected  in  the 
contours  of  mixing  ratio  in  Fig.  8b  which  have  their  greatest  vertical  extent  between  grid  points  30  and 
45. 

There  is  some  likelihood  for  the  fluid  near  and  paraUel  to  the  front  to  behave  quasi-two 
dimensional.  It  is  well  known  that  two-dimensional  turbulence  cascades  energy  toward  the  larger 
scale  (Kraichnan  and  Montgomery,  1980).  To  retain  the  true  character  of  the  flow,  three 
dimensional  turbulence  is  necessary.  In  our  study  we  have  found  diat  the  horizontal  variation  of 
meteorological  fields  in  the  vicinity  of  the  cold  front  are  very  dependent  on  the  horizontal  mixing  (7). 

In  Figs.  9a  and  9b  we  show  the  mixing  ratio  field  for  the  lowest  sigma  level,  48  hrs  into  the  forecast, 
for  simulations  made  with  and  without  the  horizontal  turbulence.  Without  the  horizontal  mixing  the 
gradients  become  too  large  and  unrealistic  (Fig.  9b).  In  our  three  dimensional  TT  formulation  we 
found  that  x=100  in  (6)  gave  reasonable  results.  The  results  in  Fig.  9b  are  seen  to  compare 
reasonably  well  with  the  control  Blackadar  simulation  shown  in  Fig.  9c.  Increasing  the  horizontal 
mixing  by  reducing  the  value  of  x  even  more  would  appear  beneficial. 

The  rms  error  at  850  mb  in  the  forecast  temperature  is  shown  in  Fig.  10.  Results  from  four 
forecasts  are  displayed.  The  control  run  is  the  Penn  State/NCAR  model  with  the  existing  Blackadar 
boundary  layer  scheme  and  K  theory  horizontal  diffusion.  The  control  was  run  with  a  Kuo  cumulus 
parameterization  and  with  an  explicit  cloud  scheme,  labeled  (clouds)  in  Fig.  10  (Hsie  and  Anthes, 
1984).  In  the  explicit  cloud  scheme  prognostic  equations  are  included  for  cloud  water  and  rain  water. 
Also,  evaporation  of  cloud  and  rain  water  in  unsaturated  layers  are  part  of  the  process  but  the  ice  phase 
is  not  considered  in  our  calculations.  In  the  TT  parameterization  version  we  run  identical  simulations 
except  no  cumulus  parameterization  was  utilized  in  any  of  our  numerical  forecasts.  Note  in  Fig.  10 
that  the  rms  errors,  determined  when  forecast  and  radiosonde  values  are  compared,  are  least  for  the 
transilient  approach  with  explicit  clouds,  except  near  the  end  of  the  72  hr  forecast.  Overall,  carrying 
the  clouds  explicitly  made  little  difference  in  these  rms  statistics,  but  we  found  that  feedback  processes, 
e.g.,  cloud  and  rainwater  evaporation,  are  of  the  utmost  importance.  Evaporation  is  included  when  the 
clouds  are  explicitly  carried;  otherwise,  the  current  version  of  the  Penn  State/NCAR  region  model  does 
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not  include  this  process.  Zhang  et  al.  (1988)  have  shown  that  including  realistic  physics,  e.g., 
evaporation  and  water  loading,  reduces  model  tendencies  for  over  development.  In  this  study  of  a 
new  turbulence  parameterization  cloud  and  rain  water  are  not  directly  included  in  the  turbulent  mixing. 
That  topic  is  left  for  future  consideration. 

The  effect  of  explicit  clouds  is  however  very  clear  in  the  mean  error  in  the  850  mb  temperatures 
shown  in  Fig.  11.  Note  that  a  change  in  sign  is  associated  with  whether  the  clouds  are  carried 
explicitly.  The  TT  scheme  gave  nearly  zero  mean  error  at  hours  12  and  24  ,  and  has  the  smallest  mean 
error  of  the  four  forecasts.  The  control  run,  with  and  without  explicit  clouds,  generally  has  much 
larger  mean  errors. 

In  Fig.  12  the  mean  sea-level  pressure  48  hrs  into  the  forecast  is  shown  for  the  control  runs,  Kuo 
cumulus  parameterization  (Fig.  12a),  explicit  cloud  (Fig.  12b),  and  for  the  TT  version  with  explicit 
cloud  (Fig.  12c)  and  for  the  verification  analysis  (Fig.  12d).  Note  that  the  low  pressure  centered  over 
the  great  lakes  is  best  represented  in  the  TT  forecast  which  has  the  correct  pressure  of  997  mbs  but 
over  a  reduced  area  as  compared  to  the  verifying  analysis  in  Fig.  12d.  The  contraction  of  the  surface 
low  in  the  TT  case  is  similar  to  the  control  case,  with  explicit  cloud  calculations,  which  has  a  central 
pressure  of  996  mb.  The  control  with  the  Kuo  scheme  has  a  990  central  pressure,  so  the  low  is 
deepened  too  much.  As  a  consequence  of  this  the  1008  pressure  contour  is  however  in  a  location 
closer  to  that  indicated  in  the  verifying  analysis.  Otherwise  over  most  of  the  remaining  area  the  control 
cases  are  slightly  better  by  a  small  margin,  but  the  presence  of  topography  is  a  complicating  factor 
making  interpretation  difficult. 

To  gain  some  idea  of  the  response  at  the  surface  we  now  turn  our  attention  to  Figs.  13a  and  13b 
showing  the  surface  slab  or  skin  temperature,  48  hrs  into  the  forecast  at  0000  UTC  24  April,  for  the 
control  case  with  explicit  cloud  calculations  and  the  equivalent  transilient  turbulence  simulation  (C3 
fixed  at  0.2).  A  comparison  of  these  two  figures  shows  that  they  are  nearly  the  same  in  the  north¬ 
eastern  portion  of  our  region  while  in  the  western  part  the  skin  temperature  differences  are  up  to  4 
degrees  in  some  locations.  We  expected  and  wanted  slightly  warmer  surface  skin  temperatures  in  the 
transilient  version.  Remembering  that  the  fluxes  are  being  reduced  by  C3  in  (13b),  which  compensates 
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for  the  molecular  layer,  means  that  the  overall  heat  flux  felt  by  the  lowest  layer  should  be  the  same 
order  of  magnitude  as  that  observed  in  the  control  case. 

Rainfall  is  strongly  correlated  with  convective  activity  and  the  vertical  velocity.  One  useful  gauge 
to  measure  development  is  the  omega  L2  norm,  i.e., 

(16) 

where  #  is  the  total  number  of  grid  points  in  the  volume.  The  variation  of  the  omega  norm  with  time 
for  the  OSCAR  IV  case  is  shown  in  Fig.  14.  The  conversion  is  that  lO'^  Pa  s''=  36  mbs  hr^.  The 
control  OT  Blackadar  scheme  shows  little  variation  with  time  while  the  TT  scheme  exhibits  significant 
changes,  i.e.,  note  the  increase  during  the  second  day  and  the  diurnal  variation  of  the  omega  norm.  A 
maximum  is  found  during  the  afternoon  hours  and  again  at  night,  (forecast  hours  0, 24, 48  and  72 
correspond  to  (XX)  UTC).  Allowing  Cj  to  vary  (0.1  <  C3  <  0.2,  labeled  Var  in  legend)  reduces  the 
norm  while  the  introduction  of  the  explicit  cloud  calculations,  including  evaporation,  also  reduces  the 
magnitude  especially  after  hour  40.  Therefore  details  in  the  surface  flux  and  viscous  layer  calculations 
can  make  a  significant  difference  in  the  fOTecast.  In  spite  of  the  enhanced  vertical  motions  rainfall  is 
only  slightly  heavier  (bias  score  increased  by  25%  for  threshold  >  2.64  cm.)  with  the  transilient 
turbulence  parameterization. 

The  cumulative  rain  and  evaporation  for  the  72  hr  forecast  is  shown  in  Fig.  1 5.  The  explicit  cloud 
simulations  are  labeled  'Exp'  while  'Reg'  refers  to  only  large  scale  rainfall  except  in  the  control  where 
the  Kuo  cumulus  parameterization  scheme  is  also  included.  Note  that  in  the  TT  forecasts  the  rainfall  is 
excessive  without  the  feedback  mechanisms  of  evaporation  and  water  loading.  Allowing  C3  to  vary, 
labeled  Var,  also  reduces  the  total  precipitation.  However  the  largest  difference  between  the  control 
and  the  new  TT  parameterization  is  found  in  the  cumulative  evaporation  columns.  This  difference  is 
primarily  due  to  the  inclusion  of  the  viscous  layer  in  the  moisture  flux  calculations  (15  ). 
Evapotranspiration  processes  have  not  been  included  in  any  of  our  simulations. 


b.  The  CAPTEX  case 

The  omega  norm  for  the  CAPTEX  case  is  illustrated  in  Fig.  16.  The  control  version  exhibits  only 
a  flat  curve,  while  using  the  TT  parameterization  (cj  =  0.2)  gives  a  large  diurnal  variation  which  peaks 
during  the  early  night  (forecast  hours  12,  36  and  60  correspond  to  0000  UTC).  It  is  very 
characteristic  of  mesoscale  convective  systems  to  exhibit  a  nocturnal  maximum  in  precipitation  (Kane 
et  al.,  1987).  This  diurnal  variability  has  also  been  noted  in  tropical  regions  by  Gray  and  Jacobson 
(1977)  and  Meisner  and  Arkin  (1987).  As  in  the  OSCAR  IV  case  the  inclusion  of  explicit  clouds 
makes  a  significant  difference.  Explicit  clouds  with  cloud  and  rain  water  evaporation  processes 
included  keep  the  three  diurnal  peaks  at  about  the  same  intensity,  while  the  maximum  grows  in  time 
when  evaporation  processes  are  not  considered.  These  findings  help  substantiate  the  conclusion 
reached  by  Zhang  et  al.  (1988),  .i.e.,  that  feedback  processes  are  important  to  keep  the  numerical 
models  from  developing  run-away  convection.  We  now  examine  this  and  other  feedback  processes 
in  greater  depth. 

It  have  been  known  for  some  time  that  many  processes,  e.g.,  radiation  and  cloud  and  rain  water 
evaporation,  are  important  and  should  be  included  to  realistically  model  the  earth's  atmosphere.  But 
when  these  processes  are  included  in  some  traditional  numerical  models  the  response  signal  is  only 
faintly  discernible  or  is  hard  to  interpret.  We  now  examine  the  role  of  feedback  mechanisms  utilizing 
our  newly  created  TT  model  with  the  CAPTEX  initial  and  boundary  data  sets.  Comparisons  are  also 
made  with  the  control  version.  A  complete  description  of  the  feedback  processes  utilized  in  the  Penn 
State/NCAR  region  model  is  found  in  Anthes  et  al  (1987).  Additional  sources  for  details  on  the 
physics  include;  Hsie  et  al.  (1984)  for  a  description  of  the  evaporative  processes,  Zhang  and  Anthes 
(1982)  present  details  about  the  surface  budget  calculations  and  Benjamin  and  Carlson  (1986)  outline 
the  long  and  short  wave  radiation  calculations.  In  Fig.  17a  we  see  that  when  cloud  and  rain  water 
evaporation  is  removed  from  the  transilient  turbulence  explicit  cloud  calculations  (TT  No  Cloud  Evap) 
the  omega  norm  is  substantially  increased  in  magnitude  as  compared  to  the  transilient  forecast  with  that 
process  included  (TT ).  With  the  removal  of  evaporation  in  the  air  convection  has  becomes  much 
more  pronounced  so  a  major  part  of  the  diurnal  variation  is  no  longer  present.  Likewise  we  see  that 


21 


I 

I 

removing  the  atmospheric  component  of  the  radiation  (TT  No  Atm  Rad)  enhances  the  omega  norm  as 
does  the  removal  of  radiative  cloud  feedback  to  the  surface  radiation  budget  (TT  No  Cloud  Rad). 

Thus,  all  these  feedback  reactions  help  control  the  degree  of  convection.  It  is  thus  clear  from  our 
numerical  results  that  the  feedback  processes  examined  here  are  very  important  to  our  72  hr  forecast. 

The  same  four  physical  feedback  tests,  also  for  explicit  cloud  calculations,  are  computed  for  the 
control  version  of  the  model.  The  omega  norm  for  each  is  shown  in  Fig.  17b.  Evaporation  has  a 
strong  signal  predominantly  during  the  latter  half  of  the  fraecast.  The  effect  of  the  other  two  physics 
mechanisms  are  more  difficult  to  interpret  since  after  72  hrs  both  omega  norms  are  less  in  magnitude 
than  the  full  physics  case  (heavy  solid  line).  Not  only  is  the  shape  of  the  curves  different  but  the 
response  and  the  timing  of  the  response  is  also  clearly  different  than  that  shown  in  Fig.  17a. 

We  postulate  that  the  feedback  or  physics  processes  become  more  important  as  the  numerical 
nradel  become  better  able  to  simulate  the  atmosphere.  The  fact  that  different  responses  are  clearly 
visible  in  Fig.  17a  within  just  twelve  hours  into  the  forecast  implies  that  numerical  weather  prediction 
and  climate  models  can  not  ignore  these  (and  other)  processes  and  expect  to  simulate  the  atmosphere 
successfully.  In  fact,  with  the  TT  parameterization  we  found  that  the  model  returns  a  coherent  signal 
in  the  omega  norm  (not  shown)  when  the  net  radiative  flux  is  changed  by  just  one  percent.  Even 
though  this  signal  is  very  faint  it  is  still  detectable  within  twelve  hours  in  the  difference  field  between 
the  forecasted  omega  norms  computed  with  100%  and  99%  of  the  net  radiative  flux. 

We  now  examine  one  of  the  most  difficult  processes  to  simulate  correctly,  i.e.,  precipitation.  The 
change  in  the  total  area  hourly  precipitation,  a  measure  of  the  precipitation  rate,  is  shown  in  Fig.  1 8  for 
the  TT  procedure  (heavy  line)  and  for  the  convective  and  large  scale  values  associated  with  the  control 
or  standard  Penn  State/NCAR  model.  Notice  the  diurnal  variation  in  the  TT  forecasted  precipitation. 
Peak  values  occur  during  the  early  night  time  hours  but  some  small  oscillations  occur  during  early 
afternoon  maximum  surface  heating  periods.  As  expected  the  precipitation  and  omega  norm  have  very 
similar  behavior  as  seen  by  comparing  Figs.  16  and  18.  In  contrast,  in  the  control  simulations  the 
Kuo  convection  scheme  peaks  at  maximum  surface  heating  times  during  the  first  half  of  the  forecast 
period  (Fig.  1 8).  The  following  36  hrs  is  a  period  of  predominantly  increasing  precipitation.  The 
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large  scale  rainfall  (control)  shows  little  variation  but  a  gradual  increase  peaks  31  hrs  into  the  forecast. 

The  question  remains;  are  there  observations  verifying  the  diurnal  early  night  time  peaks  predicted 
using  TT?  To  assist  with  answering  this  question  we  examine  satellite  visible  and  IR  images.  Four 
separate  time  periods  are  shown  in  Fig.  19.  Note  how  the  convective  cloud  activity  increases  and  then 
decays  over  the  time  period  2100,  0000, 0300  and  at  0600  UTC  for  the  25th  and  26th  of  September 
1983.  The  maximum  convective  activity  occurs  somewhere  between  between  hours  0000  and  03(X) 
UTC.  Clearly,  upward  vertical  motion  should  be  associated  with  these  centers  of  development. 

The  omega  field  predicted  in  the  TT  forecasts  are  radically  different  from  that  given  by  the  control. 
The  magnitudes  and  scale  of  the  vertical  motion  is  clearly  different  as  revealed  in  Figs.  20a  and  20b, 
representing  the  horizontal  slice  for  the  8th  sigma  level  (approximately  766  mb  standard  atmospheric 
pressure  level)  for  the  TT  and  the  control  36  hr  forecasts,  respectively.  Note  how  with  TT  (Fig.  20a) 
there  are  rolls  in  the  right  center  of  the  figure  whose  slope  parallels  the  U.  S.  East  Coast  and  the 
Appalachian  Mountains.  These  feature  are  absent  in  the  control.  Comparing  the  region  of  dashed 
lines  or  negative  omega  values,  in  Fig.  20a  with  the  cloud  cover  in  Fig.  19  gives  a  good  correlation. 
Clearly  the  magnitudes  of  these  updrafts  are  much  larger  than  those  in  Fig.  20b.  The  placement  of  the 
updrafts  in  TT  coincides  very  well  with  the  rainfall  pattern  particularly  in  Mexico  and  in  the  Colorado 
Rockies,  but  the  rainfall  amounts  associated  with  some  of  the  more  intense  updrafts,  are  much  too 
heavy. 

The  accumulated  rainfall  as  a  function  of  time  is  shown  in  Fig.  21.  Note  that  the  totals  for  the  TT 
approach  and  the  control,  with  Kuo  cumulus  parameterization,  are  nearly  the  same.  The  diurnal 
oscillation  is  apparent  ii;  the  TT  results.  In  the  control  forecast  the  convective  rainfall  is  larger  than  the 
large-scale  nonconvective  component.  Because  the  TT  approach  parameterizes  the  mixing  in  the 
vertical  from  unresolved  horizontal  scales  we  do  not  find  it  necessary  to  include  cumulus 
parameterization. 

The  variation  of  the  total  cumulative  rainfall  and  evaporation  for  Five  different  72  hr  forecasts  is 
shown  in  Fig.  22.  The  TT  forecasts  are  labeled  TT  and  those  from  the  control  are  labeled  with 
Blackad;ir.  Both  explicit  cloud  calculations  (Exp)  and  the  regular  (Reg),  which  includes  the  Kuo 
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cumulus  parameterization  for  the  control  only,  are  presented.  Note  the  significant  reduction  in  total 
rainfall  for  the  explicit  cloud  cases.  Removing  the  hcnizontal  mixing  in  the  turbulence  calculations 
(7a),  labeled  No  Horiz,  increases  the  rainfall.  Also  obvious  in  Fig.  22  is  the  significant  differences  in 
the  reccnded  amount  of  evaporation  derived  from  the  surface. 

TTie  heat  flux  for  the  72  hr  forecast  period  computed  at  point  (x2S^y2s)  is  shown  in  Fig.  23a.  (The 
numbering  of  the  grid  begins  at  the  lower  left  hand  comer.)  Values  computed  from  (13)  and  (14),  and 
those  computed  in  the  control,  are  shown.  The  control  surface  flux  is  larger  in  magnitude  but  the 
pattern  is  very  distorted  at  the  lowest  atmospheric  level  (HFLUX  15)  during  the  first  twelve  hours  of 
the  forecast,  in  contrast  to  that  given  by  the  more  physically  realistic  approached  use  in  this  study. 
Likewise  the  moisture-flux  patterns  (Fig.  23b)  are  very  similarly  behaved.  In  Fig.  23b  we  see  that 
the  control  shows  some  tendency  for  an  enhancement  of  the  moisture  flux  with  time  in  contrast  to 
nearly  cyclic  behavicM-  in  the  TT  model  calculations.  The  large  fluctuation  seen  with  the  control  are 
absent  in  our  TT  results. 

The  bias  in  the  rainfall,  measuring  the  tendency  of  the  model  to  forecast  too  small  (<1)  or  too  large 
(>1)  an  area  of  precipitation  for  a  72  hour  forecasts,  is  illustrated  in  Fig.  24.  Including  the  clouds  and 
evaporation  in  the  TT  cases  reduces  the  over-estimation  for  the  heavier  rainfall  but  does  affect  the  light 
rain  cases  because  of  the  water  loading  by  the  clouds.  Otherwise  the  TT  parameterization  with  explicit 
clouds,  even  with  the  much  larger  omega  norm,  gives  accurate  results  with  respect  to  the  placement  of 
the  precipitation,  but  the  amounts  of  precipitation  is  sometimes  overestimated  in  local  areas  by  a  factor 
of  two  or  three. 

This  accuracy  is  also  revealed  in  the  850  mb  geopotential  heights.  The  rms  error  with  re.spect  to 
radiosonde  measurements  is  shown  in  Fig.  25a.  The  forecasts  utilizing  the  transilient  turbulence 
parameterization  are  clearly  superior,  i.e.,  the  rms  error  is  reduced  by  more  than  a  third.  The  mean 
error  is  likewise  significantly  reduced  using  the  TT  scheme  as  indicated  in  Fig.  25b. 
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6.  Conclusions 

The  first-order  non-local  turbulence  parameterization  known  as  transilient  turbulence  (11)  has 
been  incorporated  into  the  Penn  State/NCAR  15  sigma  level  regional  primitive  equation  mesoscale 
model.  This  required  that  the  host  model  be  modified  by  removing  the  existing  representations  for  the 
horizontal  and  vertical  K-theory  eddy  diffusion  and  the  dry  and  moist  convective  adjustment  routines. 
In  their  place  we  substitute  the  unified  approach  described  by  the  TT  parameterization.  Additionally,  a 
sixth  order  implicit  tangent  filter  was  inserted  to  remove  numerical  noise,  thus  separating 
computational  considerations  from  the  physical  turbulence  parameterization.  In  the  surface  heat  and 
moisture  flux  calculations  a  new  procedure  incorporating  a  molecular  viscous  layer  was  used  to 
improve  the  surface  physics.  Together,  the  changes  espoused  here  represent  a  major  departure  from 
the  traditional  modelling  norm,  as  represented  by  the  host  model.  To  test  the  TT  parameterization  we 
made  72  hr  forecasts  using  the  OSCAR  IV  and  CAPTEX  initial  and  boundary  data  sets. 

In  TT  the  parameterization  is  based,  in  this  study,  on  the  turbulent  kinetic  energy  equation.  The 
vertical  component  in  the  TT  scheme  is  applied  first  and  turbulent  grid  points  are  identified.  A 
horizontal  mixing  scheme  is  then  applied  only  between  neighboring  grid  points  that  were  turbulent  in 
the  vertical.  Both  schemes  together  comprise  the  three  dimensional  TT  procedure.  The  calculations 
of  vertical  mixing  within  the  unevenly-spaced  sigma  coordinate  grids  were  modified,  as  described  in 
Section  2,  from  the  conceptual  formulation  described  in  Appendix  A.  The  gauge  used  to  measure  the 
desired  amount  of  mixing  was  precipitation.  We  found  in  both  the  OSCAR  IV  and  CAPTEX  cases 
that  the  amount  and  nature  of  the  vertical  mixing  had  a  direct  bearing  on  the  amount  and  location  of  the 
precipitation.  The  modifications  described  in  Section  2  yielded  satisfactory  results  for  our  vertical  grid 
configuration.  Additional  improvements  are  still  possible.  Nevertheless  the  current  results  from  two 
case  are  very  promising  and  show  that  major  modelling  improvements  can  result.  Additional  case 
studies  are  being  used. 

The  three  dimensional  TT  parameterization  scheme  was  found  to  be  numerically  stable,  to 
accurately  portray  boundary  layer  characteristics,  and  to  greatly  improve  the  forecast  accuracy  for  the 
CAPTEX  case.  By  far  the  largest  amount  of  turbulence  occurred  in  the  boundary  layer.  The  new 
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surface  heat  flux  calculations,  which  includes  a  viscous  layer,  improved  the  communication  between 
the  earth's  surface  and  the  surrounding  boundary  layer.  This  enabled  the  TT  scheme  to  build  the 
well-mixed  boundary  layer.  During  the  morning  and  early  afternoon  the  boundary  layer  grew  because 
of  the  surface  heating.  Thus  the  number  of  grid  points  experiencing  turbulence  reflects  the  diurnal 
cycle. 

The  omega  norm,  and  hence  the  omega  fields  themselves,  also  experience  a  diurnal  variation,  but 
now  the  maximum  occurs  during  the  early  night  time  hours.  We  believe  the  magnitude  of  this  change 
to  be  too  large  but  the  placement  of  the  updrafts  coincide  nicely  with  precipitation  verification.  The 
forecasted  precipitation  maximums  also  coincide  with  the  nocturnal  omega  maximum.  These 
nocturnal  precipitation  maximums  are  common  to  land-based  mesoscale  ccMivective  systems,  so  the 
turbulence  forecasts  are  thought  to  simulate  the  atmosphere  correctly.  Satellite  IR  images  verify  our 
forecasts  by  showing  enhanced  convective  cloud  cover  during  early  night  time  hours. 

Because  of  the  large  variation  of  the  omega  field  in  dme,  and  because  of  the  larger  omega 
magnitudes,  the  TT  model  is  much  mcoe  sensitive  to  how  physical  quantities  are  represented.  For 
example,  evaporative  cooling  now  becomes  extremely  important  Explicitly  including  clouds,  even  for 
the  80  km  horizontal  grid  considered  here,  becomes  important  because  the  evaporative  processes  are 
only  included  in  that  portion  of  the  model  code.  Atnxispheric  radiation  and  cloud  radiation  feedback 
processes  that  influence  the  surface  calculations  are  also  shown  to  be  significant  Boundary-layer 
parameters  and  calculation  details  were  also  found  to  be  important  We  suspect  that  additional 
physical  processes  like  the  role  of  evapotranspiration,  cirrus  clouds  and  cloud  micro-physics  may  now 
play  a  much  larger  role  in  the  model  with  the  turbulence  parameterization  than  in  the  control  version. 

Even  though  there  are  significant  reductions  of  error  in  the  boundary  layer  calculations  and 
improvements  in  precipitation  prediction  in  predominantly  the  CAPTEX  case,  statistically  there 
remains  little  difference  between  the  control  and  the  TT  fcaecasts  at  500  mbs  and  above.  However 
there  are  local  differences  primarily  due  to  turbulence  in  the  vicinity  of  jets,  but  for  the  geopotential 
heights  at  500  mbs  the  overall  forecasts  have  nearly  identical  rms  error  statistics,  as  determined  from 
radiosonde  reports.  This  is  true  for  both  the  OSCAR  IV  and  CAPTEX  cases. 
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In  the  Penn  State/NCAR  regional  model  the  computational  time  required  for  the  transilient 
turbulence  calculations  is  more  than  three  times  that  required  with  the  Blackadar  or  control  version  of 
the  model.  However,  with  the  new  turbulence  parameterization  the  magnitudes  in  the  forecasted 
omega  fields  are  more  than  double  that  given  by  the  control.  Because  the  horizontal  divergence  is  scale 
dependent  we  would  expect  similar  behavior  of  the  omega  magnitudes  if  the  grid  resolution  was 
increased  by  a  factor  of  two.  This  latter  calculation  would  however  take  about  eight  times  that  required 
by  the  course  grid.  Thus  the  transilient  turbulence  scheme  represents  a  desirable  alternative. 

Continued  research  and  better  model  and  algorithm  development  should  reduce  the  computational 
expense.  Clearly  the  computational  overhead  would  be  much  less,  in  a  relative  sense,  in  numerical 
models  which  use  a  larger  time  step  size,  e.g.,  semi-implicit  time  integration  or  semi-Lagrangian 
procedures.  The  TT  scheme  is  ideally  suited  for  non-staggered  grids  thus  the  Arakawa  A  grid  is  the 
ideal.  Staggered  grids  greatly  enhance  the  computational  expense  and  the  necessary  interpolations 
unduly  modify  the  presence  and  concentration  of  the  turbulence. 

We  have  shown  that  the  transilient  turbulence  parameterization  improves  the  boundary  layer 
behavior,  predicts  CAT  and  forecasts  omega  fields  which  exhibit  a  large  diurnal  variation  with  a 
nocturnal  maximum.  Forecast  statistics  are  definitely  improved  for  the  CAPTEX  case.  We  also 
found  that  physical  feedback  processes  like  atmospheric  radiation,  cloud  radiative  properties  on  the 
surface  calculations  and  cloud  and  rainwater  evaporation  all  give  an  important  signal  with  a  positive 
contribution  to  the  72  hr  forecasts.  These  processes  are  shown  to  be  all  important  in  simulating  the 
atmosphere  accurately. 

We  believe  that  with  the  aid  of  additional  refinements  the  transilient  turbulence  procedure  has  the 
potential  to  significantly  improve  numerical  weather  prediction. 
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Appendix  A 

Transilient  Parameterization  for  Unequal  Grid  Spacing 

a.  A  physical  interpretation  of  nonlocal  mixing 

First,  picture  a  column  of  10  equally-spaced  grid  points  as  sketched  in  Fig  Ala.  Let  each  grid  box 
be  filled  with  one  unit  mass  of  air.  Furthermore,  assume  for  the  sake  of  example  that  there  is  total 
exchange  of  only  those  units  masses  of  air  as  sketched  in  Fig  Ala.  For  example,  all  of  the  air  from 
grid  box  a  moves  into  box  d,  and  all  of  the  air  originally  in  d  moves  to  a,  with  similar 
interpretations  for  the  other  mixing  arrows  drawn.  The  transilient  matrix  for  this  case  is  symmetric,  as 
shown  in  Fig  Alb. 

Next,  group  those  unit  boxes  of  air  into  larger  grid  boxes  of  different  sizes  (indexed  as  1  to  4),  as 
shown  in  Fig  Ale,  but  assume  the  same  mixing  processes  are  occurring.  The  new  grid  box  1  is  very 
small,  and  after  mixing  it  holds  precisely  one  unit  mass  that  came  from  the  new  grid  box  3;  thus  C13  = 

1 .0.  Grid  box  3  is  larger,  thus  the  unit  mass  of  air  that  returned  from  box  1  accounts  for  only  33%  of 
the  final  mass  in  box  3;  that  is,  C31  =  0.33.  Similar  arguments  can  be  used  to  determine  the  other 
transilient  coefficients,  which  arc  sketched  in  Fig  Aid.  Thus,  the  transilient  matrix  for  unequally- 
spaced  grid  points  is  asymmetric,  even  when  symmetric-type  exchange  mixing  processes  are 
occurring. 

The  definition  of  a  transilient  coefficient  is  a  very  specific  "destination"-oriented  definition:  Cjj  is 
the  fraction  of  air  in  destination-grid-box  i  that  came  from  source-grid-box  j.  This  definition  is  still 
valid  for  unequally-spaced  grid  points  and  must  be  interpreted  strictly;  namely,  Cy  is  NOT  the  fraction 
of  air  in  source-box  j  that  left  to  go  to  box  i.  The  first  index  always  represents  the  destination,  the 
second  always  represents  the  source,  and  the  value  of  c  is  always  taken  with  respect  to  the  destination. 
Such  a  "destination"  definition  must  be  used  in  order  for  the  prognostic  equation  to  work: 

Si(t+At)  =  Zj  Cij(At)  Sj  ,  which  weights  the  various  incoming  tracer  values  by  their  relative 
contributions  to  the  final  total  mass  in  the  destination  box. 
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With  the  strict  definition  of  the  transilient  coefficients,  each  row  of  the  transilient  matrix  must  sum 


to  unity: 


1 


(Al) 


Physically,  this  is  interpreted  as  air-mass  conservation  in  the  destination  grid  boxes.  After 
mixing,  each  grid  box  is  still  filled  with  the  same  mass  of  air  that  it  started  with,  and  this  air  had  to 
have  come  from  somewhere.  Fig  Aid  exemplifies  this  interpretation. 

However,  each  column  (which  represents  each  source)  does  not  necessarily  sum  to  unity  because 
of  the  "destination"  definition  of  the  transilient  coefficients.  If  mj  and  mj  represent  the  total  masses  of 
air  within  grid  boxes  j  and  i,  respectively,  then  the  fraction  of  air  originally  in  source-box  j  that  left  to 
go  to  box  i  is  (mj/mj)  cy.  Conservation  of  source  material  (ie.,  state  conservation)  is  now: 
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m. 

m7 
i  =  1  J 


(A2) 


as  is  exhibited  in  Fig  Aid.  Both  (A  1)  and  (A2)  are  valid  for  any  mixing  process,  even  if  there  is  not 
an  equal  exchange  between  the  upper  and  lower  grid  points. 

For  the  special  case  of  perfect  mass  exchange  as  shown  in  Fig  Al,  the  coefficients  in  the  upper 
triangle  of  the  matrix  are  related  to  those  in  the  lower  triangle  by; 


m. 

c..  =  — -  c.. 
J>  m.  u 
J 


(A3) 


where  the  first  index  of  c  always  represents  the  destination.  The  mixing  potential  parameterization 
used  in  the  body  of  this  paper  does  not  discriminate  between  upward  and  downward  mixing,  and  thus 
is  described  by  (A3). 
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b.  A  parameterization 

The  approach  taken  above  to  change  from  a  large  number  of  equal-size  grid  boxes  to  a  smaller 
number  of  unequal-size  boxes  can  be  also  used  to  parameterize  the  transilient  coefficients  as  a  function 
of  mixing  potentials,  Yjj.  Consider,  for  example,  the  potential  for  mixing  between  each  of  the  unit 
masses  within  box  4  (source)  and  the  one  unit  mass  within  box  1  (destination).  The  potential  for 
mixing  between  the  unit  mass  g  and  unit  mass  a  is  Ygg  (=Yga).  Similarly,  the  potentials  between  the 
other  unit  masses  in  4  and  the  unit  mass  in  1  are  Y^j, ,  Ygp ,  and  Ygq  .  The  total  mixing  potential 
between  box  I  and  all  the  unit  masses  in  box  4  is  the  sum  of  all  the  individual  potentials:  Ygg  +  Y^h  + 

Yap  +  Yaq. 

Mixing  potentials,  however,  are  determined  by  wind  and  temperature  differences  between  the  grid 
boxes  according  to  (1).  For  the  case  of  the  column  of  10  smaller  grid  boxes,  we  expect  that  each  of 
the  individual  muting  potentials  ( Yah  .  etc.)  could  differ  from  each  other  corresponding  to  different 
winds  and  temperatures  in  those  regions.  However,  for  the  case  of  the  unequally-spaced  coarser  grid 
box,  we  know  only  one  temperature  and  wind  for  the  whole  box,  because  anything  smaller  is 
unresolved  subgrid  scale.  Thus,  Ygg  =  Yah  =  Yap  =  Yaq  =  Y 14.  As  a  result,  the  total  mixing  potential 
between  boxes  1  and  4  is:  Ygg  +  Yah  ■*"  Yap  +  Ygq  =  4  Y14.  In  general,  the  total  mixing  potential  is 
mj  Yjj,  where  we  see  that  the  mass  of  the  source  box  is  a  weighting  factor  for  the  mixing  potential. 

A  similar  weighting  does  NOT  occur  for  the  destination  boxes.  For  a  coarse  grid  box,  such  as 
destination-box  4,  100%  of  the  air  mixing  into  any  one  of  the  smaller  destination  unit-mass  boxes 
contributes  only  25%  to  the  total  mass  in  box  4.  Thus,  the  total  contribution  to  the  overall  mixing  into 
box  4  is  0.25  Yga  +  0.25  Yha  +  0.25  Ypa  0.25  Yqg  .  Because  each  of  the  mixing  potentials  for  the 
smaller  unit  masses  contained  within  box  4  are  equal  to  Y41  (i.e.,  no  subgrid  differences),  we  find 
that  the  potential  for  mixing  into  box  4  is  only  Y41,  with  no  weighting  by  the  mass  of  the  destination 
box. 

The  row  norm,  RNj,  defines  the  total  mixing  potential  within  a  row  [see  (2)  in  Section  2],  and 
must  include  the  source-mass  weights  as  shown  above.  The  single  maximum  row  norm  defines  the 
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maximum  total  mixing  potential  that  can  occur  within  a  contiguous  patch  of  turbulence,  and  is  our 
definition  of  a  matrix  norm,  HYII^  [see  (3)  in  Section  2]. 

Our  basic  assumption  for  turbulence  closure  is;  that  the  fraction  of  air  mixed  into  a  destination  box 
from  some  source  is  proportional  to  the  fraction  of  total  mixing  potential  associated  with  that 
source/destination  pair  [described  by  (4a)  in  Section  2].  Because  the  row  norms  for  each  row  (except 
one)  are  smaller  than  the  matrix  norm,  that  implies  that  those  rows  (destinations)  have  less  total 
turbulent  mixing  into  them  than  the  one  lucky  destination.  As  a  result,  the  no-mixing  transilient 
coefficient  (i.e.,  cy)  is  larger  in  those  rows  such  that  the  sum  of  each  row  of  transilient  coefficients  still 
equals  unity  [see  (4b)  in  Section  2]. 
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APPENDIX  B 


Implicit  Tangent  Filter 

The  set  of  symmetric  low-pass  implicit  tangent  filters  of  order  2p  are  defined  by 

[S2P]U„P  -I-  (-l)Pe[L2p]u^F=  [S2p]u^  ,  (Bl) 

Here  u„^  is  the  filter  variable  while  is  the  unfiltered  quantity.  The  parameter  e  is  the  filter 
factor.  In  (B 1)  L2unis  the  finite-difference  analog  of  the  second  derivative  of  u  at  grid  location  n 
multiplied  by  the  square  of  the  grid  step  size,  i.e.,  L2un  =  u„.]  -  2Un  +  u„+i.  Similarly,  order 
2p  of  the  operator  is  analogous  in  the  same  way  to  a  derivative  of  order  2p.  The  coefficients 
associated  with  the  (L)2p  operation  are  identical  to  those  in  the  binomial  expansion  of  (a  -  b)2p . 
Also  S2un=!  Un.i  +  2u„  +  Un+.i ;  Consequently  (S/2)2  may  be  thought  of  as  an  averaging  operator. 
The  coefficients  for  (S)2p  are  the  coefficients  in  the  row  with  2p+l  entries  in  Pascal's  triangle. 
The  amplitude  response  of  the  implicit  tangent  filter  of  order  2p  is  given  by 

F(k)  =  [\+e  tan2p(K6/2)]-i .  (B2 ) 

Here  5  is  the  grid  step  distance  while  K  is  the  wave  number. 


39 


Figure  Captions 


Fig  1 .  Schematic  of  the  vertical  grid  spacing  of  the  sigma  and  half-sigma  levels  in  the  host  Penn 

State/  NCAR  mesoscale  model.  Temperature,  moisture,  and  horizontal  winds  are  at  the  half¬ 
sigma  levels,  while  vertical  velocity  arc  at  the  full  sigma  levels. 

Fig  2.  Schematic  (a)  highlighting  one  element  of  a  transilient  matrix,  and  (b)  indicating  the 

associated  nonlocal  mixing  between  grid  boxes  in  a  column.  The  vertical  distance  between 
grid  points  2  and  4  is  also  indicated. 

Fig  3.  Schematic  vertical  cross-section  through  a  host  riKxlel  demonstrating  vertical  and  horizontal 
turbulent  mixing.  The  subset  of  grid  boxes  that  participate  in  vertical  mixing  arc  indicated 
with  (o).  Those  points  that  arc  also  horizontal  neighbors  with  other  turbulent  grid  points  can 
participate  in  horizontal  mixing  (shaded). 

Fig  4.  Schematic  horizontal  cross-section  through  a  host  model  showing  the  horizontal-mixing 
stencil.  For  any  grid  point  in  the  model  (solid  circle  labeled:  o),  it  is  allowed  to  mix  with 
only  the  subset  of  four  or  less  immediate  neighbors  that  are  turbulent  (solid  circles  labeled:  1 
to  4). 

Fig.  5.  The  filter  response  verses  wave  length  nAx  after  720  applications  of  the  implicit  tangent  filter 
and  the  4th  order  K  theory  diffusion  in  the  Penn-State/NCAR  regional  model.  The  shaded 
region  indicates  the  difference  in  response  when  using  the  minimum  or  maximum  value 
allowed  for  K. 

Fig.  6.  The  percent  of  the  grid  points  turbulent  as  a  function  of  time  and  pressure. 
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Fig.  7.  The  horizontal  distribution  of  turbulence  at  (a)  ct=.94  and  (b)  at  o=.74  twenty  one  hours  into 
the  forecast. 

Fig.  8.  Vertical  cross  sections  of  the  (a)  potential  temperature  (K)  and  (b)  mixing  ratio  (kg  kg'^. 
Turbulence  indicated  by  stipple. 

Fig.  9.  Mixing  ratio  (kg  kg‘1)  for  the  15th  sigma  half  level  at  forecast  hour  48  for  (a)  TT  with 

horizontal  mixing,  (b)  TT  without  horizontal  mixing  and  (c)  for  the  control  Penn  State/NCAR 
model  forecast.  Contour  interval  is  0.1  E-2  while  the  labels  are  scaled  by  l.OE+4. 

Fig.  10.  The  rms  error  in  the  850  mb  temperature  (K)  for  the  control  forecast,  for  the  standard  Penn 
State/NCAR  package  of  the  Blackadar  boundary  layer  formulation  with  K  theory,  with  the 
Kuo  cumulus  parameterization  and  with  explicit  clouds.  Also  illustrated  is  the  TT  results 
with  and  without  clouds. 

Fig.  11.  Mean  error  for  the  cases  displayed  in  Fig.  7. 

Fig.  12.  Mean  sea  level  pressure  (mb)  for  (a)  the  control  run  with  Kuo  cumulus  parameterization  and 

(b)  with  explicit  clouds.  The  TT  parameterization  scheme  with  explicit  clouds  is  shown  in 

(c) .  The  verifying  48  h  analysis  is  given  in  (d).  Contour  interval  is  4  mb. 

Fig.  13.  The  forecasted  skin  temperature,  at  hour  48,  for  the  control  (a)  and  for  (b)  the  TT 
formulation.  Both  are  for  explicit  cloud  calculations.  Contour  interval  is  3  K. 

Fig.  14.  Omega  norm  as  a  function  of  time.  Note  the  difference  between  the  control  and  the  TT 

forecasts.  Influence  of  allowing  the  viscous  layer  coefficient  C3  to  vary  in  time  (Var)  and  of 
explicit  cloud  calculations  (clouds)  are  also  illustrated. 
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Fig.  15  Comalative  rain  and  evaporation  fw  the  72  hr  forecasts  with  (Exp)  and  without  (Reg)  explicit 
clouds.  Explicit  cloud  calculations  ccmtain  the  influence  of  the  feedback  mechanisms  of 
rainwater  evaporation  and  water  loading.  Allowing  C3  to  vary  (V ar)  reduces  the  rainfall  in 
the  TT  fcnecasts. 

Fig.  16  Omega  norm  as  a  function  of  time  for  the  CAPTEX  case.  Note  the  diurnal  variation  for  the 
TT  parameterization  cases. 

Fig.  17  Omega  norm  (CAPTEX)  showing  the  influences  of  feedback  processes.  In  (a)  the  effects  on 
the  TT  forecasts  are  illustrated  while  (b)  outlines  the  results  in  the  control  forecasts. 

Fig.  18  The  change  in  the  hourly  forecasted  precipitation  is  shown  for  the  TT  forecast  (CAPTEX) 
and  fcH*  the  convective  and  large  scale  components  of  precipitation  in  the  control  forecast. 

Fig.  19  Satellite  IR  images  showing  cloud  cover  at  (a)  2100  UTC  25  September  1983  and  for  (b) 
0000 ,  (c)  0300  and  (d)  0600  UTC  26  September  1983. 

Fig.  20  Horizontal  fields  of  omega  (CAPTEX)  at  the  8th  sigma  level  36  hrs  into  the  forecast  ((X)00 
UTC  26  September  1983)  for  (a)  TT  and  (b)  for  the  control.  Contour  interval  is  0.15E-3 
while  labels  are  scaled  by  0. 1 E+6.  Dashed  lines  indicate  upward  nx>tion. 

Fig.  21  The  accumulated  rainfall  is  shown  as  a  function  of  time  for  the  TT  forecast  and  the  large  scale 
and  COTvcctive  components  of  the  precipitation  in  the  control  forecast. 

Fig.  22  Cumulative  rainfall  and  evaporation  for  TT  with  (Exp)  and  without  explicit  clouds  (Reg), 
with  no  horizontal  mixing  (No  Horiz),  and  for  the  control  with  Kuo  cumulus 
parameterization  and  with  explicit  clouds. 
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Fig.  23  Comparison  of  the  heat-flux  (a)  and  moisture-flux  (b),  given  in  units  of  W  m-2  between  TT 
and  the  control  at  grid  location  (X25,y25)-  TTie  TT  heat  and  moisture  flux  is  for  the  surface 
and  sigma  level  15  while  for  the  control  we  present  both  values. 

Fig.  24  Bias  scores  in  rainfall  (CAPTEX)  measuring  the  accuracy  of  the  coverage  of  the  rainfall  after 
72  hours. 

Fig.  25.  The  850  mb  (a)  rms  and  (b)  mean  error  in  the  geopotential  heights  for  the  control 
(Blackadar)  and  TT  approach,  with  and  without  explicit  clouds. 

Fig  A1 .  (a)  Schematic  of  a  column  of  10  equally-spaced  grid  points  labeled  a  -  q,  each  representing  a 
unit  mass  of  air  contained  within  the  grid  boxes  drawn.  For  sake  of  demonstration,  assume 
that  air  is  completely  removed  from  some  grid  boxes  and  replaced  by  air  from  other  grid 
boxes,  as  indicated  by  the  arrows. 

(b)  Symmetric  transilient  matrix  corresponding  to  (a).  Note  that  this  matrix  is  shown  upside- 
down,  so  that  each  row  of  the  matrix  corresponds  to  the  associated  height  of  the  grid  box  of 
part  (a).  The  main  diagonal  is  highlighted  by  a  shaded  line,  and  zero  elements  are  left  blank. 
Each  row  and  column  of  the  matrix  sums  to  one. 

(c)  Schematic  of  a  column  of  4  unequally-spaced  grid  points  based  on  groupings  of  the  unit 
masses  from  (a).  The  same  amount  of  mass  exchange  as  (a)  is  assumed  to  occur,  as 
indicated  by  the  arrows. 

(d)  Asymmetric  transilient  matrix  corresponding  to  (c).  Each  row  sums  to  one,  but  each 
column  does  not. 
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FIGURE  16 
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